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Abstract 

We discuss and illustrate through numerical examples the relations between generating 
functionals, thermodynamic consistency (in particular the virial-free energy one), and 
uniqueness of the solution, in the integral equation theory of liquids. We propose a new 
approach for deriving closures automatically satisfying such characteristics. Results from 
a first exploration of this program are presented and discussed. 
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I Introduction 



Integral equation theories (lET) of the hquid state statistical mechanics are valuable tools for 
studying structural and thermodynamic properties of pairwise interacting fluid systems [1,2]. 
Many of these approximations to the exact relation between pair potential and pair correlation 
functions have been proposed in the last half century, starting from the pioneering works [3-5] 
to the most refined and modern approximations [6-10] which may approach the accuracy of 
computer simulation with a negligible computational cost. 

The functional method in statistical mechanics [1] provides the most general and sound 
starting point to introduce lET as approximations of the exact functional relations and it is 
the classical statistical mechanics counterpart of the quantum density functional theory. 

Notwithstanding the success of present lET to describe the structure of simple one com- 
ponent systems, considerable work is still devoted to derive improved approximations which 
could accurately describe the thermodynamics as well. Also applications to non simple or 
multicomponent systems are still subject of current studies. 

Actually, the description of thermodynamics is one weak point of lET approaches: rea- 
sonable and apparently harmless approximations to the potential-correlation relations usually 
result in a dramatically inconsistent thermodynamics where many, if not all, among the exact 
sum rules derived from statistical mechanics, are violated. 

The problem of thermodynamic inconsistency, i.e. the inequivalence between different routes 
to thermodynamics, actually plagues the lET approach to the point that the degree of incon- 
sistency between different formulae for the same quantity is used as an intrinsic measurement 
of the quality of a closure. 

In the past, some discussion of the thermodynamic consistency appeared in the literature. 
Hyperncttcd chain approximation (HNC) was recognized as a closure directly derivable from an 
approximation for the free energy functional [11] , thus exhibiting consistency between the virial 
formula and the thermodynamic expression for the pressure. However, this limited consistency 
is not enough to guarantee a unique and faithful description of the phase diagram. Apart the 
problem of the remaining inconsistencies, the descriptions of the critical points and spinodal 
lines arc seriously inadequate. 

Extensive work on HNC [12-14] showed that in place of a true spinodal line, it is more 
appropriate to describe the numerical results as due to a region in the thermodynamic plane 
where no real solution of the integral equation exists. In particular, Belloni[12] showed that 
the disappearance of the solution originates from a branching point where two solutions merge, 
instead than from a line of diverging compressibility. Thus, we have direct evidence that HNC 
may have multiple solutions, at least in part of the phase diagram. 

Empirical improvements on HNC have been proposed [6, 9, 10] providing in many cases ex- 
cellent results for one-component simple fluids. However, although reduced, the thermodynamic 
inconsistency problem remains and the multiple solution problem is completely untouched. 

In this work we start an investigation of a new approach to lET directly addressing the 
two points of uniqueness of the solution and thermodynamic consistency. The basic idea is 
to constrain the search for new closures within the class of generating functionals which are 
strictly convex free energy functionals, thus enforcing the virial-energy consistency as well as 
the uniqueness of the solution. 

In particular, in the present paper we try to answer the following questions: i) does at 
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least one strictly convex free-energy functional of the pair correlation function exist? ii) what 
is the nature of the resulting spinodal line (if any), iii) what is the quality of the resulting 
thermodynamic and structural results? iv) does the simultaneous requirement of consistency 
and uniqueness automatically provide improved results? 

As we will show, we have a positive answer for i) , a thorough and interesting characterization 
for ii), some interesting indications for iii), and a partly negative answer for iv). 

However, we can show that it is possible to exploit the control provided by the generating 
functional approach to easily generate new closures and we feel our procedure could be the 
basis of a more systematic approach to lET. 

In section II we recall the connections between closures, generating functionals, thermody- 
namic consistency and uniqueness of solutions and we illustrate them in the well known case 
of HNC approximation. In section 111 we introduce two straightforward extensions of HNC 
intended to cure its problems. In Section IV numerical results are presented and discussed. In 
section V we show two possible improvements of the closures studied. 



II Thermodynamic consistency and uniqueness of the 
solution of integral equations 

Since the work by Ohvares and McQuarrie [15] it is known the general method to obtain 
the generating functional whose extremum with respect to variations of the direct (c(r)) or 
total {h{r)) correlation functions results in the closure relation, provided the Ornstein-Zernike 
equation is satisfied. 

For example, if we have a closure of the form 

p'c{r)^nKr),mr)} , (2-1) 
where 0(r) is the pair interaction potential and ^ is an arbitrary function, the functional 

- J drh{r) dt'if{th{r),P(f){r)} + constant^ , (2.2) 



is such that the extremum condition 

SQ 



5h{r) 

is equivalent to 



= , (2.3) 



p'^h{r)^^{h{r),P(f){r)} + p J h{\r-r'\)^{h{r'),P(f){r')}dr' . (2.4) 

Olivares and McQuarrie also showed how to find the generating functional if the closure is 
expressed in the form 

p'^h{r) = ^{c{r),/3(j){r)} . (2.5) 
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In appendix A we discuss the extension of their method to the case of a closure written as 



p2c(r) = ^{7(r),/?0(r)} , (2.6) 

where 7(r) = h{r) — c(r) is the indirect correlation function. Notice that most of the modern 
closures correspond to this last case. 

The possibility of translating the original integral equation into an extremum problem allows 
to get an easy control on two important characteristics of the approximation: thermodynamic 
consistency between energy and virial routes to the thermodynamics and uniqueness of the 
solution. 

Indeed, once we get the generating functional Q, due to the approximations induced by the 
closure, there is no guarantee that its value at the extremum is an excess free energy. In order 
to be a free energy, the functional should satisfy the condition 

-,9{r) , (2.7) 



50(r) 2' 

where g{r) = h{r) + 1 is the pair distribution function. 

Even if this condition is not new, and mention to it is present in the literature [16], we 
discuss it in appendix B as well as its consequences on the thermodynamic consistency between 
the virial pressure and the density derivative of the free energy. 

Another issue where the generating functional approach is useful is the problem of multiple 
solutions of the integral equations [12]. In particular, the analysis of the convexity properties 
of the generating functional is a very powerful tool [17, 18]. 

Let us illustrate this techniques in the case of HNC closure. It is well known [11, 15] that 
the HNC equation with closure 

c(r) = h{r) - In [g{r)e^'^^''^ , (2.8) 

can be derived from the variational principle 



5h{r) 



= , (2.9) 



where 



with 



^Fihj^J'ozihj+THNcih] , (2.10) 



^^{ph{k)-Hi + ph{k)]} , 

^ J'HNcih] ^p'Jdv{l+ g{r) [in (^(r)e^'^W) - l] - /i^(r)/2} 



(2.11) 



Let us call h(r) the extremum of JF, solution of the variational principle (2.9). It can be 
shown (see appendix B) that, within an additive constant, J-'[h]/{2(3p) is the excess Helmholtz 
free energy per particle of the liquid. This ensure thermodynamic consistency between the 
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route to the pressure going through the partial derivative of the free energy and the one going 
through the virial theorem (see appendix B). In addition, it allows to get a closed expression 
for the excess chemical potential without further approximations [19, 20]. This feature is highly 
desirable for applications of lET to the determination of the phase diagrams. 

Moreover if we can prove that T ^ defined on some convex set of trial correlation functions 
is a strictly convex functional, then we know that if a solution to (2.9) exists, it corresponds 
to a minimum and is unique. A functional T is strictly convex if for all |/(r) e and |/(r) ^ 0, 
we have 



/ 



We calculate the second functional derivatives as follows 



6h{r)6h{r') ' J (27r)3 [l + ph{k)Y ' 

(2.13) 



- r') hxx - 1 



5^:FHNc[h] _ 

5h{r)5h{r') "'^ ^ ' \g{r) 
Recalling that the static structure factor S{k) = 1 + ph{k), we find for A 

A/P'- [dry\r)('-l] . (2.14) 



(27r)3 52(A;) J ^ ' ' \g{r 

Now, the most interesting results would be to show the strict convexity of the HNC functional 
over the convex set of all the admissible pair correlation functions (all the h{r) > —1 and 
properly decaying to zero at large distance. 

However, this is not the case for HNC. It has not been possible to show the positive definite- 
ness of equation (2.14) and it has been shown [12] that in some region of the thermodynamic 
plane HNC does exhibit multiple solutions. 

The best we can do is to obtain a more limited result. Calling gi = sup g{r) {gi > 1 is the 
height of the first peak of the pair distribution function) and using Parseval theorem, we find 

from which we deduce that A > on the following set of functions 

D=\^h{r)\Q<S{k) < y/g,/{g,-l) VA;} . (2.16) 

We conclude that JF defined on any convex set of functions Dc G D is strictly convex. Near the 
triple point we arc sure we arc out from such set since the first peak of the pair distribution 
function for the Lennard- Jones fluid is ~ 3 [21], so that \/'gi/{gi'--T) ~ 1.2. The first 
peak of the static structure factor is also close to 3. Then we are not inside D and the HNC 
approximation may have multiple solutions [12]. 

Instead, if we are in the weak coupling regime, the previous conditions tells us that there is 
a range where the branch of solutions going to the perfect gas limit is unique and quite isolated 
from other solutions. 
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Ill Extensions of HNC 



The generating functional approach can be used in a systematic way to look for better closures. 
We think that this way, we can obtain a less empirical search method for improving closures. 

In the following we report some preliminary analysis we have done. As a first test of our 
program, wc have restricted our investigations to simple modifications of HNC functional. As 
we will discuss later, such a choice is certainly not optimal. However, we can learn enough to 
consider the approach worthwhile of further investigations and we feel the results are interesting 
in order to reveal more details about the characteristics of the solutions of the highly non linear 
lET. 



111.1 The HNC/H2 approximation 

We want to modify the HNC closure in order to have an integral equation with a generating 
functional which is strictly convex without having to restrict its definition domain. We choose 
as our modified HNC (HNC/H2) closure [22] 

c(r) = h{r) - Hg{r)] - (3(t){r) - ah\r) , (3.1) 

with a a parameter to be determined. The new closure generating functional is 

y'HNC/H2[h] =p'Jdr{l + g{r) [in (^(r)e^^M) - l] - h\T)/2 + ah^r)/^] . (3.2) 

Its second functional derivative with respect to h is 

(3.3) 

Recalling that h — g — 1 and g{r) > for all r, we see that for a = 1/2 

--l + 2ah= ^ ^ > . (3.4) 

9 9 

Then Thnc/H2 is a convex functional and since Toz is unchanged and strictly convex (see 
appendix C), their sum, the generating functional of the integral equation, is strictly convex. 

Moreover {J^oz[h] + ^HNC/H2[h]}/{'^Pp) continues to be the excess Helmholtz free energy 
per particle of the liquid since equation (2.7) holds (see appendix B). 

We have then an integral equation which is both thermodynamically consistent (the pressure 
calculated from the virial theorem coincides with that one calculated from the Helmholtz free 
energy) and with a solution which, when it exists, is unique. 

111.2 The HNC/H3 approximation 

In the same spirit as in subsection lll.l we can try to add a term in the HNC/H2 closure 
c(r) = h(r) - ln[g{r)] - (3(t){r) - ah^{r) - -fh^(r) , (3.5) 



S J^HNC/H2[h] 

Sh(r)Sh(r') 



= p'5(r - r' 



1 + 2ah{r) 
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with a and 7 parameters to be determined. We call this approximation HNC/H3. The closure 
generating functional is 

Thnc/hM =p'Jdr{l + g{r) [in (^7(r)e^<^M) - l] - h\r)/2 

+ah^{r) /3 + 7/i^(r) /4} . (3.6) 
Its second functional derivative with respect to h is 

1 



Sh{r)5h{r') ^ ^ ' 



- 1 + 2ah{r) + 3-fh\r) 



= p^S{r-v')^-^{l-2agir) + 3^gir)[l-gir)]} . (3.7) 

In order to have the right hand side of this expression positive for g > the only choice we 
have is to set q; = 1/2. In this way 

(l-^)[l-2a^ + 37^(1-^)] = (1-^)^(1 + 37^) , (3.8) 

and we see that J-'hnc/hs is a convex functional if we additionally choose 7 > — l/[3sup5'(r)]. 

Once again {J^oz[h] + J^HNC/H3[h]}/{'^Pp) is the excess Helmholtz free energy per particle 
of the liquid and the thermodynamic consistency virial-free energy is ensured. 



IV Numerical results 

To solve numerically the OZ plus closure system of nonlinear equations we used Zerah' s algo- 
rithm [23] and Fourier transforms were done using fast Fourier transform. In the code we always 
work with adimensional thermodynamic variables T* — l/(/3e), p* — pa^, and P* — Pa^/e, 
where a and e are the characteristic length and characteristic energy of the system respectively. 
We always used 1024 grid points and a step size Ar = 0.025o". 

The thermodynamic quantities were calculated according to the statistical mechanics for- 
mulae for: the excess internal energy per particle 

poo 

C/^^7Ar = 27rp / (j){r)g{rydr , (4.1) 
Jo 

the excess virial pressure 

(3P^Ip-1^-It,(3p j^'^girydr , (4.2) 
the bulk modulus calculated from the compressibility equation 

PXt S{k = 0) 

where xt is the isothermal compressibility, and the bulk modulus calculated from the virial 
equation 

Qpv 



dp 
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For the calculation of Bp once g{r) and c(r) had been calculated, Lado' s scheme for Fourier 
transforms [24] was used to determine dg{k)/dp. Even if slow, this allows us to explicitly 
calculate and later invert the coefficients matrix of the linear system of equations which enters 
the calculation of dg{k)/dp. 

IV. 1 Inverse power potentials 

The general form of the inverse power potential is 

where 3 < n < oo. For this class of fluids the thermodynamics depends only from the dimen- 
sionless coupling parameter 

^ = (p(7VV2)(/3e)^/" . (4.5) 

We performed our calculations on the n = 12, 6, 4 fluids at the freezing point. We compared 
three kind of closures: the one of Rogers and Young [25] (RY) with thermodynamic consistency 
virial-compressibility and known to be very close to the simulation results, the hypernetted- 
chain (HNC) closure, and the HNC/H2 described in subsection III.l. In each case we compared 
our data with the Monte Carlo (MC) results of Hansen and Schiff [26]. 

IV. 1.1 The inverse 12th power potential 

The freezing point for this fluid is at z = 0.813. The RY a parameter to achieve thermodynamic 
consistency at this value of z is 0.603. 

In table i we compare various thermodynamic quantities (the excess internal energy per 
particle, the excess virial pressure, the bulk moduli) obtained from the RY, the HNC, and 
the HNC/H2 closures. In the MC calculation the excess virial pressure is 18.7 and the bulk 
modulus 72.7. 

In flgure 1 we compare the MC, the HNC, and the HNC/H2 results for the pair distribution 
function. 

IV. 1.2 The inverse 6th power potential 

The freezing point for this fluid is at z = 1.54. The RY a parameter to achieve thermodynamic 
consistency at this value of z is 1.209. 

In table ii we compare various thermodynamic quantities (the excess internal energy per 
particle, the excess virial pressure, the bulk moduli) obtained from the RY, the HNC, and 
the HNC/H2 closures. In the MC calculation the excess virial pressure is 38.8 and the bulk 
modulus 110.1. 

IV. 1.3 The inverse 4th power potential 

The freezing point for this fluid is at z = 3.92. The RY a parameter to achieve thermodynamic 
consistency at this value of z is 1.794. 
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In table iii we compare various thermodynamic quantities (the excess internal energy per 
particle, the excess virial pressure, the bulk moduli) obtained from the RY, the HNC, and 
the HNC/H2 closures. In the MC calculation the excess virial pressure is 108.7 and the bulk 
modulus 156. 

In figure 2 we compare the MC, the HNC, and the HNC/H2 results for the pair distribution 
function. 

IV. 2 The spinodal line 

In this subsection we study a pair potential with a minimum In particular we chose the Lennard- 
Jones potential 



where e and a are positive parameters. The critical point for this fluid is at [27] T* = 1.3120 ± 
0.0007, p* = 0.316 ± 0.001, and PI = 0.1279 ± 0.0006. 

Integral equations usually fail to have a solution at low temperature and intermediate densi- 
ties, i.e. in the two- phases unstable region of the phase diagram. In particular it is well known 
that the HNC approximation is unable to reproduce the spinodal line, the locus of points of 
infinite compressibility in the phase diagram [12]. This is due to the loss of solution as one 
approaches the unstable region on an isotherm from high or from low densities. The line of 
loss of solution, in the phase diagram, is called termination line. According to the discussion 
of section II, the loss of solution for the HNC approximation can be traced back to the loss 
of strict convexity of the generating functional [28]. Indeed, using HNC approximation, we 
computed the bulk modulus from the compressibility equation Be, on several isotherms as a 
function of the density. At low temperatures we found that both at high density and at low 
density we were unable to continue the isotherm at low values of Be- Zerah' s algorithm either 
could not get to convergence or it would converge at a non physical solution (with a pole in the 
structure factor at some finite wavevector k). Since HNC/H2 has, by construction, an always 
strictly convex generating functional, we expect it to be able to approximate a spinodal line 
(there should be no termination line). 

In Figure 3 we show the behavior of Be on several isotherms as a function of density, 
calculated with the HNC/H2 approximation. We see that now there are no termination points. 
Be never becomes exactly zero and the low temperature isotherms develop a bump in the 
intermediate density region. The same plot for the bulk modulus calculated from the virial 
pressure Bp, shows that at low temperatures this bulk modulus indeed becomes zero along the 
isotherms both at high and low densities. 

In figure 4 the pressure is plotted as a function of the density on several isotherms for the 
HNC/H2 approximation. Apart from the fact that we find negative pressures, the isotherms 
have a van der Waals like behavior. 

The graphical analysis of the pressure plotted as a function of the chemical potential shows 
that the coexistence of the two phases (points where the curve crosses itself) is possible and is 
lost between T* — 1.1 and T* — 1.2. There generally are two points of coexistence. 




(4.6) 
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V Improving the closures 



The numerical results for HNC/H2 exhibit interesting features as far as the coexistence re- 
gion is concerned but show unambiguously a worst agreement with the MC structural data in 
correspondence with a marginal improvement in the thermodynamics. 

We feel that the main problem is the difficulty of an accurate description of the bridge 
functions in terms of powers of the pair correlation function. Recent investigations on improved 
closures seem to point to the indirect correlation function 7(r) or some renormalized version of 
it, as the best starting point for progress. However, before moving to more complex relations or 
functional dependences, we have explored two possible directions for improving the HNC/H2 
closure. In the first approach we have tried to follow the MHNC approach by Lado et al. 
[29] . In the second we have explored the possibilities of optimization offered by the numerical 
coefficient of the cubic term in the generating functional. 

V.l Pseudo-bridge functions for HNC/H2 

From the graphical analysis of the pair distribution function it is known [1] that g{r) may be 
written as 

g{r)^eM-(3(f>{r) + j{r) + B{r)) , (5.1) 

where 7(r) = h{r) — c(r) is the sum of all the series type diagrams and B{r) the sum of bridge 
type diagrams. If we take 

B{r)^-^h\r) + G{r) , (5.2) 

we have that our HNC/H2 approximation amounts to setting G{r) = 0. Roscnfeld and Ashcroft 
[6] proposed that B{r) should be essentially the same for all potentials 0(r). We now make a 
similar proposal for the G function and we will refer to it as pseudo bridge function. In the 
same spirit of the RHNC approximation of Lado [29] we will approximate G{r) with the G 
function of a short range (reference) potential (t)o{r). Assuming known the properties of the 
reference system, we can calculate the G function as follows 

G'o(r) = In [^o(r)e^^°«] - 7o(r) + \hl{r) . (5.3) 

The reference HNC/H2 (RHNC/H2) approximation is then 

^(r)=exp(-/30(r) + 7(r)-^/i2(r) + Go(r)) . (5.4) 

An expression for the free energy functional can be obtained turning on the potential (j){r) 
in two stages: first, from the noninteracting state to the reference potential 0o('") and then 
from there to the full potential 0(r). To this end we write 

0(r;Ao,Ai) = Ao0o(r) + AiA0(r) , (5.5) 
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with A0(r) = 0(r) — 0o('^)- Following the same steps as in [30] we obtain for the excess free 
energy per particle 

r"^ = /i + /2 + /f +A/3 (5.6) 
where the first two terms were already encountered in section II 

Pfi = \pjdv{l-r g{r) [In (^(r)e^'^W) - l] - h\r)/2 + h\r)/&] , (5.7) 

The third term is assumed known 

= -lpJd^l' dXMr; Ac, 0) ^'^'^^'j = /3(/^°^ - /f ^ - /f ) , (5.9) 

here /(°^ is the excess free energy per particle of the reference system and f[^\ /2°^ are defined as 
in equations (5.7), (5.8) for the reference potential and its corresponding correlation functions. 
The last term is 

/3A/3 = -^p J dr J^' dXMr; 1, ^O^^^J;^ ■ (5-10) 

According to our proposal, G is insensitive to a change in potential from 0o to 0. We may then 
approximate this last term as follows 

(3Afs ^ ~p j drGo{r)[g{r) - ^o(r)] . (5.11) 

Now that we have the free energy we may consider it as a functional of both h{r) and Goir) 
and take its variation with respect to these functions. We find, 

l^^jexc ^ ^pj dr{c{r)-h{r)+h\r)/2 + ]n[g{r)e'^'t'^''^] - Go{r)} Sh{r) 

-\p j dr[g{T) - g^{r)]5Go{r) . (5.12) 

It follows that the free energy is minimized when both the RHNC/H2 closure (equation (5.4)) 
is satisfied and when the following constraint 

j dr[g{T)-go{T)]5Go{r)^Q , (5.13) 

is fulfilled. 

Taking the second functional derivative of /^^"^ with respect to h{r) we find that also this 
free energy is a strictly convex functional of the total correlation function. This property was 
lacking in the RHNC theory and constitutes the main feature of the RHNC/H2 closure. As 
already stressed in section III.l it ensures that if a solution to the integral equation exists it 
has to be unique. 
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The constraint, as for RHNC, gives a certain thermodynamic consistency to the theory (see 
[30]). If we choose a hard sphere reference potential (f)o{r) = (j)o{r;a) which depends on the 
length scale a, the optimum values of the parameters that makes the generating functional a 
free energy can be determined by the constraint (5.13) which becomes 



However, neither the hard-sphere pseudo bridge functions nor some empirical attempt to 
model the unknown function via a Yukawa function provided useful results. 

V.2 Optimized HNC/H3 approximation 

For 7 = HNC/H3 reduces to HNC/H2. For 7 > the first peak of the pair distribution 
function is dumped respect to the one of the pair distribution function calculated with HNC/H2. 
For 7 < the first peak increases giving in general a better fit to the simulation data. 

In figure 5 we compare the pair distribution function of the Lennard- Jones fluid near its 
triple point, calculated with a molecular dynamic simulation [21], the HNC/H2 approximation, 
the approximation HNC/H3 with 7 = —0.203 (at lower values of 7 Zerah' s algorithm would fail 
to converge), and the approximation HNC/H3 with 7 = —0.1 (when the generating functional 
of HNC/H3 is still strictly convex). As we can see HNC/H3 fits the simulation data better 
than HNC/H2 even if the first peak is still shghtly displaced to the left of the simulation data, 
a well known problem of the HNC approximation [6] . 

The best results are given by HNC/H3 with 7 = —0.203. Note that the HNC/H3 generating 
functional at this value of 7 is not strictly convex (strict convexity is lost for 7 < —1/9). The 
first peak of the static structure factor is at ka ~ 6.75 and has a magnitude of 2.41, a quite low 
value for a liquid near the triple point. We have calculated the pressure and the internal energy. 
We found /3P/p ~ 3.87 and C/^^7(A^e) ~ -5.72 (very close to the HNC results pP/p ~ 3.12 
and U^^'^/{Ne) ~ —5.87) to be compared with the simulation results [31] 0.36 and —6.12 
respectively. The bulk moduli are — 11.74 and Bp ~ 36.61 which shows that at the chosen 
value of 7 we do not have the thermodynamic consistency virial- compressibility and we do not 
improve on HNC inconsistency (using HNC we find Be ~ 7.09 and Bp ~ 32.72). 

VI Conclusions 

In this paper we have analyzed the relations between generating functionals, thermodynamic 
consistency and uniqueness of the solution of the integral equations of liquid state theory. We 
think that the requirement of deriving from a free energy and the uniqueness of the solution are 
two important ingredients to enforce in the quest for better closures. The former requirement 
is of course crucial to get virial-energy consistency. But it is also important to get integral 
equations able to provide a closed formula for the chemical potential without additional ap- 
proximations. This last issue looks highly desirable for applications of lET to the determination 
of phase diagrams. The latter is certainly a useful constraint from the numerical point of view 
but it is also a very strong condition, probably able to avoid some non physical behavior in the 
coexistence region, although this point would deserve further investigation. 




(5.14) 
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In this work, we have started an exploration of the capabihties of the combined requirement 
of consistency and uniqueness, starting with simple modifications to the HNC closure, corre- 
sponding to the addition of a square and a cubic power of h{r) in the HNC functional. We 
found a couple of approximations (HNC/H2 and HNC/H3), which have built in the virial-free 
energy thermodynamic consistency and have a unique solution. 

We numerically tested these closures on inverse power and the Lennard- Jones fluid. Prom 
the tests on the inverse power potential fluids one can see that the HNC/H2 approximation 
is comparable to HNC for the thermodynamic quantities and performs worst than RY and 
even HNC for structural properties. The tests on the Lennard-Jones fluid revealed as this 
approximation does not suffer from the presence of a termination line (present in HNC and 
almost all the existing closures). This allowed us to follow isotherms from the low density 
to the high density region and this behavior would be very useful in the study of the phase 
coexistence. However, the thermodynamic results show only a marginal improvement on HNC 
and the structure is deflnitely worse. 

Our trials to improve HNC/H2 in the same spirit of the modifled HNC approaches did not 
succeed. We feel that the main reason is in the difficulty of modeling the real bridge functions 
through a polynomial in the function h{r). In this respect, approaches based on generating 
functionals depending on the indirect correlation function 7(r) look more promising but we 
have not tried them yet. 

Much better results for the structure are found with HNC/H3 as is shown in figure 5. 
However, probably for the same reasons just discussed, one has to renounce to have an approx- 
imation with a strictly convex generating functional depending on h{r). The thermodynamics 
reproduced by HNC/H3 is not yet satisfactory: due to the slight left shift of the main peak 
of the g{r) the calculated pressure misses the simulation result. Nonetheless the presence of 
the free parameter 7 in HNC/H3 leaves open the possibility of imposing the thermodynamic 
consistency virial-comprcssibility. If the value of the parameter needed to have the consis- 
tency is bigger than — 1/[3 sup ^'(r)] then we would have an approximation which is completely 
thermodynamically consistent and have a unique solution. This strategy may eventually lead 
to discover that the price we have to pay to have a completely thermodynamically consistent 
approximation is the loss of strict convexity of the generating functional. 

A Appendix: Generating functionals of 7 

Often in the numerical solution of the OZ plus closure integral equation use is made of the 
auxiliary function 7(r) = h{r) — c(r). Suppose that the closure relation can be written as 



where ^ is a local function of the function 7 and has a dependence on the value of the pair 
potential not explicitly shown. 

We want to translate the integral equation into a variational principle involving functionals 
of 7(r). Then we introduce a closure functional J^dil] such that 



P'c(r) = -*{7(r)} 



(A.1) 



*{7(r)} 



(A.2) 
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and an OZ functional J-'oz,c[y] such that, when c(r) and 7(r) satisfy the OZ equation, we have 



(A.3) 



Then when both the closure and the OZ relations are satisfied, the functional T — ^oz,c 
is stationary with respect to variations of 7(r), i.e. 



5m 



= 



(A.4) 



This is the variational principle sought. 

Now, we want to find Toz,c- The OZ equation in k space is 

pc^{k)^ p^{k)c{k)-^{k) = ^ . 

When we solve it for c we find two solutions 

. -f ± Vf 2 + 4f 



2p 



(A.5) 



(A.6) 



where f (A;) = P7(^) is always positive since 



r = p^hc = 



(A.7) 



5'(A;) being the liquid static structure factor which is positive definite for all k. Since c{k) is 
a function which oscillates around 0, where c is negative we have to choose the solution with 
the minus sign, where it is positive the one with the plus sign. In particular, if the isothermal 
compressibility of the liquid xt, is smaller than the one of the ideal gas Xtj have that 



c(0) 



1 - 



XT 



< 



(A.8) 



and we have to start with the minus sign. 

The functional we are looking for is then (see equation (30) in [15], with the constant set 
equal to zero) 



^OZ,cb] 



dt / dr7(r) 



(27r)3 2 



tT{k) + Sc{k)Jt^T^{k) + 4tT{k) 



(A.9) 



where Sc{k) is +1 when c{k) > and —1 when c{k) < 0. Rearranging the integrals and making 
the change of variable y — tV we find 



-P/A + s,{k) 



I 



(l + f/2) ^(l + f/2)'-l 
-In (^l + r/2 + ^J(l + r/2y -1^ I . (A.IO) 
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If the closure relation has the form 

p^h{r) = -^{7(r)} , (A.ll) 

we can derive the corresponding functional using the same procedure. The final result is a 
functional J-oz,h['y] which differs from (A. 10) for a plus sign in front of the first term in the 
integral. 

However, by examining their second functional derivatives, we notice that both J^oz,c[y] 
and J-oz,h[l'] ai'c not certainly convex or concave. Thus, any check of the convexity properties 
of generating functionals of the j{r) function should be done on the full functional. 



B Appendix: Thermodynamic consistency 



For a homogeneous liquid interacting through a pair potential 0(r), the Helmholtz free energy 
per particle / can be considered a functional of (p. Indeed, in the canonical ensemble, one has 



= /^/o - ^ In ( ^ / exp -P^- 0( 



dvi • • - dr 



N 



:b.ii 



where /o is the free energy per particle of the ideal gas (0 = 0) and V is the volume of the 
liquid. Taking the functional derivative with respect to (5(f){r) one finds 



air) 



(B.2) 



Imagine that we found a functional .4.([/i], [0],p, /?) that has an extremum for those corre- 
lation functions that solve the OZ and the closure system of equations. Suppose further that 
such functional has the following property 



SpA _ p 
5p(t){r) ~ 2 

which can be rewritten more explicitly as follows 



5I3A 



9{r) , 



(B.3) 



(B.4) 



Evaluating this expression on the correlation function h solution of the OZ plus closure system 
of equations, which is an extremum for A, we find 



S(3A 



5(5(t){r) 



(B.5) 



Then we can write 



M([;i],[0],p,/3)= / dv 



6PA 



P<l>{r)+V{[h],p,P) , 



(B.6) 



[h],P,P 
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with V a functional independent of 0. Changing variables to adimensional ones, r = v*p 
and using equation (B.5) we find 

/5^([/i*],[0],p,/5) = i J dr*-g*{r*)mr*p-'/')+V{[h*],p,P) , (B.7) 

where we defined new distribution functions g*{r*) = g{r*p^^/^). If T> has no exphcit depen- 
dence on p then one readily finds 

p 9mmMp,(3) ^ fdv*nr'W{r*p-^'y*p-^i^ 

dp 6 J 

= /3P"^7p , (B.8) 

where again we used the fact that A has an extremum for h = h. We used a prime to denote 
a derivative with respect to the argument and P^^'^ is the excess pressure of the liquid. 
If T> has no explicit dependence on /3 we also find 



dpAi[h*],[cP],p,f3) 



^ 2 / ^^9{r)(f>{r) 



op 

= C/^^VAT , (B.9) 

where U'^^'^ is the excess internal energy. 

If V has no exphcit dependence on both p and p, V{[h*], p, P) = V{[h*]), we conclude from 
equations (B.8) and (B.9) that 

^([r],[0],p,/3) = /-'=(p,/?) + constant , (B.IO) 

where /^^'^ is the excess free energy per particle of the fiuid. Under these circumstances we 
see from equation (B.8) that we have thermodynamic consistency between the route to the 
pressure going through the partial derivative of the free energy and the route to the pressure 
going through the virial theorem. 

C Appendix: Strict convexity of J-'ozlh] 

It can be proven that the functional 

/r/k - 
^^{ph{k)-Hi + ph{k)]} , (c.i) 

defined on the convex set 

D^^{h{r)\S{k)>0 yk} , (C.2) 

is a strictly convex functional. The strict convexity is a trivial consequence of the strict con- 
vexity of the integrand in equation (C.I). 
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It remains to prove that Dc is a convex set. Given two elements of this set h' and h" , we 
need to show that h — Xh' + {1 — X)h" is an element of D(. for all A e [0, 1]. Since 

S{k) = l + ph{k) 

= 1 + p[Xh'{k) + (1 - X)h"{k)] 

= l + X[S'ik)-l] + il-X)[S"ik)-l] 

= XS'{k) + (1 - X)S"{k) > VA e [0, 1] , (C.3) 

then Dc is a convex set. 
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closure 




pp{^}/p - 1 


Be 


Bp 


RY {a = 0.603) 


2.626 


18.359 


69.782 


70.125 


HNC 


3.009 


21.036 


45.278 


80.430 


HNC/H2 


3.200 


22.372 


52.661 


87.255 



Table i: We compare various thermodynamic quantities as obtained from the RY, the HNC, 
and the HNC/H2 closure, for the inverse 12th-power fluid at the freezing point {z = 0.813). 
jjexc I jVg) is the excess internal energy per particle, (iP^""^ j p - 1 the excess virial pressure. 
Be and are the bulk moduli calculated from the compressibility and the virial equations 
respectively. 
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closure 




pp(^)/p-l 


Be 


B, 


RY {a = 1.209) 


4.114 


39.027 


110.952 


111.420 


HNC 


4.235 


40.178 


84.016 


113.733 


HNC/H2 


4.283 


40.635 


88.289 


115.757 



Table ii: We compare various thermodynamic quantities as obtained from the RY, the HNC, 
and the HNC/H2 closure, for the inverse 6th-power fluid at the freezing point {z = 1.54). 
jjexc I jVg) is the excess internal energy per particle, iiP^""^ j p - 1 the excess virial pressure. 
Be and are the bulk moduli calculated from the compressibility and the virial equations 
respectively. 
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closure 




pp(^)/p-l 


Be 


Bp 


RY {a = 1.794) 


8.001 


104.664 


250.106 


242.948 


HNC 


8.047 


105.277 


223.328 


244.212 


HNC/H2 


8.068 


105.542 


226.966 


257.678 



Table iii: We compare various thermodynamic quantities as obtained from the RY, the HNC 
and the HNC/H2 closure, for the inverse 4th-power fluid at the freezing point {z = 3.92). 
jjexc I jVg) is the excess internal energy per particle, iiP^""^ j p - 1 the excess virial pressure. 
Be and B^ are the bulk moduli calculated from the compressibility and the virial equations 
respectively. 
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Figure 5: R. Fantoni and G. Pastore 
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